-- methods for the simluation in the phase space

-- laws of motion
  local function newton(positions)
    local accelerations = 0*positions
    for a=1, #positions do
      for b = a+1, #positions do
        local diff = positions[a]-positions[b]
        local norm = math.sqrt(Vector.dot(diff,diff))
        if norm<1 then
          diff = 0*diff
          norm = 1
        end
        norm = norm*norm*norm
        accelerations[a] = accelerations[a] - (Masses[b]/norm) * diff
        accelerations[b] = accelerations[b] + (Masses[a]/norm) * diff
      end
    end
    return accelerations
  end
  local function deriv(state)
    return Vector{state[2], newton(state[1])}
  end
  

  -- simple euler integration
  local function euler(dt, state, deriv)
    return state + deriv(state)*dt
  end

  -- 4th order runge kutta
  local function rk4(dt, state, deriv)
    local k1 = deriv(state)
    local k2 = deriv(state + k1*(dt/2))
    local k3 = deriv(state + k2*(dt/2))
    local k4 = deriv(state + k3*(dt))
    local k = (dt/6) * (k1 + k2*2 + k3*2 + k4)
    return state + k
  end

  -- runge-kutta-fehlberg
  local function rkf(dt, state, deriv)
    local k1 = deriv(state)
    local k2 = deriv(state + dt*(k1/4))
    local k3 = deriv(state + dt*(3/32*k1 + 9/32*k2))
    local k4 = deriv(state + dt*(1932/2197*k1 -7200/2197*k2 + 7296/2197*k3))
    local k5 = deriv(state + dt*(439/216*k1 -8*k2 +3680/513*k3 -845/4104*k4))
    local k6 = deriv(state + dt*(-8/27*k1 + 2*k2 -3544/2565*k3 + 1859/4104*k4 -11/40*k5))
    local k  = dt*(25/216*k1 + 1408/2565*k3 + 2197/4104*k4 -1/5*k5)
    local kk = dt*(16/135*k1 + 6656/12825*k3 + 28561/56430*k4 -9/50*k5 +2/55*k6)
    return state +k, state+kk
  end
minscale = 100
maxscale = 0
  -- adaptive rkf
  -- not working yet, do not use
  local function rkfa(dt, st, deriv, eps)
    local s, control = rkf(dt, st, deriv)
    local err, worst = s-control, 0
	
	for _, v in ipairs(err[1]) do
		local magnitude = Vector.dot(v,v)
		if magnitude > worst then worst = magnitude end
	end

	for _, v in ipairs(err[2]) do
		local magnitude = Vector.dot(v,v)
		if magnitude > worst then worst = magnitude end
	end
	
    local scale = math.sqrt(math.sqrt(dt*eps/worst/2))
    if minscale>scale then minscale = scale end
    if maxscale<scale then maxscale = scale end
    return s
  end

  function advance(state, delta)
    return rkfa(delta, state, deriv,0.01)
	--local _,s = rkfa(delta, state, deriv,1)
    --return s
  end
  
  function advance_simple(state, delta)
    return euler(delta, state, deriv)
  end

